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The Minkowski operators (addition and substraction of sets in vectorial spaces) has been exten- 
sively used for Computer Graphics and Image Processing to represent complex shapes. Here we 
propose to apply those mathematical concepts to extend the Molecular Dynamics (MD) Methods 
for simulations with complex-shaped particles. A new concept of Voronoi- Minkowski diagrams is 
introduced to generate random packings of complex-shaped particles with tunable particle round- 
ness. By extending the classical concept of Verlet list we achieve numerical efficiencies that do not 
grow quadratically with the body number of sides. Simulations of dissipative granular materials 
under shear demonstrate that the method complies with the first law of thermodynamics for energy 
balance. 

PACS numbers: 45.70.-n 45. 40.-f 47.11. Mn 



I. INTRODUCTION 

It has passed more that 50 years since Loup Verlet pub- 
lished the milestone paper which gives birth to the Molec- 
ular Dynamics (MD) method [1]. MD is the paradigm of 
modelling complex systems as a collection of particles 
interacting each other. The method has been success- 
fully applied for computer simulations in different areas 
of Physics, Chemistry and Biology. Even engineers have 
also recognized the potential of this method, after the 
pioneer work of Peter Cundall in the extension of MD 
to model dissipative particle materials. Cundall's ideas 
on the use of discrete approach for modeling rocks and 
soils were proposed in the appendix of his PhD thesis in 
1971 This appendix led to a new area of numerical 
analysis of engineering problems, which is today known 
as Discrete Element Method (DEM). Most of the ad- 
vances in the areas of MD and DEM have been based on 
the development of efficient and robust methods to ac- 
count particle interaction. The advances in modeling of 
the real particle shape, however, has been rather slower, 
in part because the existing methods require too much 
computational effort, and partly because the difficulty to 
achieve physical correctness in the interactions. Indeed 
most of the current models use agglomerates of spheres 
to represent particles with complex shape 0, 0] . 

On the other hand. There is another fast growing area 
of research in computer graphics and rendering. New de- 
velopments are fuelled by computer games industry and 
special effects companies. There is therefore a natural in- 
terest to put these development to a good use of solving 
scientific problems. The aim of this paper is to merge the 
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molecular dynamics approach with some key elements of 
computer graphics: The Minkowski operators, Voronoi 
diagrams and Euclidian distance formulas. This combi- 
nation gives rise to a new molecular dynamics method 
that accounts both particle shape and interactions, and 
at the same time keep a reasonable balance between ac- 
curacy and efficiency. 

Contrary from the computer games, scientific and en- 
gineering applications of numerical simulations require 
high precision and interaction models which comply with 
the conservation principles of physics. Simulations with 
large number of particles has been implemented using 
three different methodologies: The first one corresponds 
to the so-called event-driven methods, where the interac- 
tion of the particles is handling via collisions [5]. These 
methods are suitable for loose particle materials, but it 
cannot be applied for dense systems because they cannot 
handle permanent or lasting contacts The second 
methodology is the Contact Dynamic methods In 
this method the equilibrium equations of the system is 
solved in each time step. The method is suitable for rest- 
ing contacts with infinite stiffness, but simulations are 
computationally expensive and lead in some cases to in- 
determinacy in the solution of contact forces [8]. This 
indeterminacy is removed by the method of Soft Molecu- 
lar Dynamics [9] . In this method the bodies are allowed 
to interpenetrate each other and the force is calculated in 
terms of their overlap. The method has been successfully 
applied for circular and spherical particles. However, the 
determination of contact force for non-spherical particles 
is still heuristic, computational inefficient, and in some 
cases it lacks physical correctness |10|, [ll| . 

Soft Molecular Dynamics Methods has been proposed 
for two-dimensional (2D) simulations using polygons 
prl , [iH. The interaction between the polygons is cal- 
culated in terms of the overlap area. The main drawback 
of this approach is that it does not provides energy bal- 
ance equations, as the elastic force between the polygons 
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does not belong from a potential The approach 

of deriving the interaction in terms of the overlap area 
is are also extremely difficult to extend to 3D, because 
the calculation of the overlap between two polyhedra is 
computationally very expensive. A solution of this diffi- 
culties has been proposed recently using spheropolygons 
flS^ . These are figures that are obtained using the math- 
ematical concept of Minkowski sum, first introduced by 
Liebling and Pourning in the modeling of granular mate- 
rials with non-spherical particles p^ . 

This paper presents an extension of the original ideas 
of Verlet, Cundall, Pournin and Liebling to model partic- 
ulate materials with complex particle shape. In Sectionllll 
we begin introducing kew elements of Computer Graph- 
ics. Then we introduce the new method of Voronoi- 
Minkowski diagrams to generate packings of complex- 
shaped particles. Section IIIII extends the concept of 
Verlet list for neighbor identification between the par- 
ticles. In section IIVI we perform numerical simulations 
of sheared granular materials with elastic, frictional con- 
tacts. We deal with energy balance and numerical effi- 
ciency. In Section [Vl we present the perspectives of this 
work to extend the model to 3D simulations. 



II. COMPUTATIONAL GEOMETRY 

Here we introduce some mathematical concepts which 
are useful to represent particle shape. After define the 
Minkowski operators and the Voronoi diagrams, we apply 
these two concepts to generate random packings of non- 
spherical particles. 

A. Minkowski operators 

Minkowski addition and substraction of sets belong to 
the area of Mathematical morphology. This is a tool 
for image processing originally developed for quantitative 
description of geological data [15]. The primary applica- 
tion of morphology occurs in by applying expanding and 
shrinking operations on binary or gray level images. Here 
we present the definition of these operations in Euclidean 
spaces and it application to the generation of objects with 
complex shape with tunable particle roundness. 



deal with the spheropolygons, which are Minkowski sum 
of a polygon with a disk, see Fig. [TJ Other exam- 
ples of Minkowski sum operations are the spherocyllinder 
(sphere line segment) |16], the Minkowski cow (non- 
convex polygon disk) [17], the spherosimplex (sphere 
simplex) UM and the spheropolyhedron (sphere poly- 
hedron) [18]. 




FIG. 1: Dilation of a square by a disk. 



2. Erosion 

The erosion of a set A by the set B is defined over an 
Euclidean space E by: 

AeB = {x\BsCA}, (2) 
where is the translation of B by the vector i.e., 

B^ = {y^x\yeB}. (3) 

Then the erosion of A by B can be understood as the 
locus of points reached by the center of B when B moves 
inside A. If A is a polygon and B is a disk of radius r, the 
erosion is a polygon inside A whose borders lie a distance 
r from A, see Fig. [2l 




FIG. 2: Erosion of a square by a disk 



i. Dilation 

Given two sets of points A and B in an Euclidean space, 
their Minkowski sum, or dilation, is given by 

A® B = {x^y \ X e A, yeB} (1) 

This operation is geometrically equivalent to the 
sweeping of one set around the profile of the other with- 
out changing the relative orientation. This paper will 



3. Opening 

The opening of A by B is obtained by the erosion of A 
by B, followed by dilation of the resulting image by B: 

AoB = {AqB)®B. (4) 

In the case of a polygon, the opening by a disk is the 
polygon with rounded corners, see Fig. [31 The degree of 
roundness increases as the radius of the disk increases. 
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FIG. 3: Opening of a square by a disk. 



B. Calculation of the mass properties of 
spheropolygons. 

Here we introduce an analytical method to calculate 
the mass, center of mass and moment of inertia of the 
spheropolygons. Consider the general case of Fig. [H 
the general spheropolygon is divided in some constituting 
parts. 




FIG. 4: General spheropolygon divided in special regions for 
the calculation of its mass properties. 



The division is as follows, first in the center (white re- 
gion in Fig. [4j) we have the original polygon that will 
be dilated by the disk element. Then we have in black 
some rectangular sections of thickness equal to the dila- 
tion radius. Finally the disk sectors in gray have an angle 
equal to f3 = 6 — ir with 6 the internal angle between the 
corresponding two segments of the original polygon. 

The spheropolygon area, and therefore its mass, is eas- 
ily calculated as the addition of the area of each sector. 
The area of rectangles is simply Ar = I x R where / is 
the length of the side and R the spheroradius. The area 
of the circular sectors is Ac 
central polygon is given by 



§R^. And the area of the 



n— 1 
i=0 



(5) 



with n the number of sides and = (x^, yi) the coor- 
dinates of the i-th vertex. 

To calculate the center of mass of the spheropolygon we 
start calculating the center of mass of each sector. The 
rectangular sectors are the easiest for this task, the center 
of mass is in the middle of them. The circular sectors 
have a center of mass in the bisector line a distance dcm = 



^^^^3/3^^^^ from the center of the circle. The coordinates 
for the central polygon center of mass are 



1 

6A 



n-l 

E 
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Xi+iVi) (6) 



-j^ n— 1 



+ yi+i){xiyi^i - Xi^iyi). (7) 



i=0 



The total center of mass is the weighted sum of all these 
center of mass coordinates. 

Now we calculate the moment of inertia. Since this 
is a 2D problem, the moment of inertia is just a scalar 
quantity. First we calculate the moment of inertia around 
the center of mass of all the regions. For the rectangles it 
is known that Ir = j2Mr{P -\- R^) with the rectangle 
mass. The moment of inertia for the circular sector is 



Mci?^(16(cos(/3)-l)+9/3^) 



(8) 



Finally for the polygon moment of inertia the formula is 



(9) 



The total moment of inertia is given by the addition of 
all these different quantities with the parallel axis theo- 
rem. This theorem states that the moment of inertia of a 
body with respect to an axis separated a distance L from 
its center of mass is 



(10) 



where Icm is the moment of inertia of the body referred 
to its center of mass and M is its mass. With this the 
mass properties calculations are complete. 



C. Voronoi-Minkowski diagrams 

Here we detail the method to generate random pack- 
ings of spheropolygons with a tuning void ratio and parti- 
cle roundness. This approach uses the concept of Voronoi 
diagrams. This is a special decomposition of the Eu- 
clidean space by space filling polygons. These polygons 
are generated by choosing a set of points in the space, 
which are called the Voronoi sites. Each site p has a 
Voronoi cell, which consists of all points closer to p than 
to any other site. 

The cells of the Voronoi construction are polygons in 
contact with each other without void spaces. Cells with 
void spaces and roundness will be constructed by using 
the concept of Voronoi-Minkowski diagrams. They are 
collections of spheropolygons resulting from the applica- 
tion of Minkowski operations to each Voronoi Cell. The 
operations are chosen such that each spheropolygon lies 
entirely in the voronoi cells, so that the spheropolygons 
has no overlap. 
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The implementation of the Voronoi-Minkowski cehs 
in 2D consists of two steps. To begin we generate the 
voronoi cehs. We set the Voronoi sites randomly dis- 
tributed in a plane. Then the polygons are constructed 
by searching the set of closest points to a particular 
Voronoi point than to any other. The resulting Voronoi 
cells has a random distribution of size and shape, as 
shown the Fig. [5l 




FIG. 5: A typical voronoi construction. The Voronoi sites 
are randomly distributed in the plane and the polygons are 
the set of closest point to each sites. 

The second step generates an spheropolygon Wi in 
each Voronoi cell Ai. We first erode the cell by a disk Bd 
of radius d. Then we apply the morphological dilation 
on the eroded polygon by a disk Br. 

W, = {AQBd)®Br. (11) 

The condition r < d guarantees that the spheropoly- 
gon is a subset of the Voronoi cell. If r = the Eq. (pTj) 
reduces to the opening defined above. In this case the 
spheropolygons touch with the neighbors by their bor- 
ders. 

The implementation of the erosion of each Voronoi cell 
is described as follows: For each vertex of the cell we 
find a sub-polygon on the inside of the voronoi polygon. 
Each vertex Xe of this sub-polygon satisfies the criterion 
that each one of them is exactly at a distance d of two 
different sides of the voronoi polygon. As can be seen in 
Fig. [6] the point Xe is obtained as the vector: 



Xe =x^dcot{0/2)e2 ^ d{k x 62). (12) 

Where x is the position of the vertex of the polygon, 62 
is the unitary vector of the edge E2, is the angle between 
the edges Ei and E2 and /c is a vector perpendicular to 
the plane. Since the set of polygon vertices are arrange in 
a counterclockwise orientation the product k x 62 points 
always to the inside of the polygon. As shown the Fig. 
[71 the set of points satisfying Eq. ([12]) gives us a possible 
vertex of the erosion sub polygon. After we chose only 
those points which are not closer than d to any other 
polygon side, as shown in Fig. [71 




FIG. 6: The erosion of the voronoi polygon is done by finding 
the sub polygon vertices considering that they are at distance 
d of two sides of the polygon. Ei and E2 are the vector 
representing two different edges of the voronoi polygon. The 
dashed vector goes from the intersection to the vertices of the 
sub polygon 




FIG. 7: The initial voronoi polygon A (left) is eroded by 
a disk element Bd of radius d. The eroded polygon A Q Bd 
(center) is then dilated by the same disk element Bd producing 
the desired spheropolygon {A Bd) Bd (right). Note the 
pathological case of points 1 and 2, since they are closer than 
d to one of the polygon sides, they are substituted by the 
point 3. 

Once the operation given by Eq. ([TT]) is done for all 
the voronoi cells, we can construct Voronoi-Minkowski 
diagrams like the one showed in Fig [HI The parameter 
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d of Eq. (pT]) controls the initial packing of the system. 
The parameter r is a measure of the angularity of the 
spheropolygons: as r decreases the particles became more 
angular. Therefore the Voronoi-Minkowski diagrams al- 
low us to generate random packings of particles with tun- 
able void ratio and roundness. 



middle point of the overlap region between the vertex 
and the edge: 




FIG. 8: An array of voronoi spheropolygons constructed by 
the method described above. The disk element used for the 
dilation process is slightly smaller than the one used in the 
erosion in order to allow void space between particles. In the 
opening the two disk elements are the same and there are no 
void spaces between sides shared by two polygons. 



III. INTERACTION FORCE 

To solve the interaction between spheropolygons we 
consider all vertex-edge distances between the polygons 
base. Let's take two spheropolygons Wi = Pi ^ Br^ and 
Wj = Pj ®Brj . We will call Pi and Pj the polygons base, 
and ri and rj the sphero-radii. Each polygon is defined 
by the set of vertices {Vi} and edges {Ej}. The force Fij 
acting on the spheropolygon i by the spheropolygon j is 
defined by: 



4 = -FJ^ = E W,^.) + E ^(^.-'^0, (13) 

where F{V, E) is force between the vertex V and the 
edge E. The torque resulted from this force is given by 



= j:v.EAmuEj)-n)xF{Ei,V,) 



(14) 



R{V, E)=X + {ri + h{V, E)) 5 \ 
2 ' '\\Y-X\ 



where 6{V,E) is defined as 



(15) 



SiV,E) = {n + rj-d{V,E)), 



(16) 



and d{X^E) = HY" — X|| is the distance from the vertex 

V to the segment E. Here X is the position of the vertex 

V and Y is its closest point on the edge E. The ramp 
function (x) returns x if x > and zero otherwise. 

The calculation of the force is extremely simple. It does 
not need neither calculations of overlap areas between nor 
calculation of Minkowski operators. Further reduction of 
the number of calculation is also possible by restricted the 
force calculation only on those vertex-edge pairs which 
are relative close each other. 

With this aim we introduce the neighbor list^ which is 
the collection of pair particles whose distance between 
them is less than 2a. The parameter a is equivalent 
to the Verlet distance proposed by Verlet to speed up 
simulations with spherical particles . 

A link cell algotithm [ll| is used to allow rapid calcu- 
lation of this neighbor list. First we decompose the space 
domain of the simulation into a collection of square cells 
of side D -\- where D is the maximal diameter of the 
particles. Each cell has a list of particles hosted in it. 
Then the neighbors for each particle are searched only in 
the cell occupied by this particle, and its eight neighbor 
cells. 

For each element of the neighbor list we create contact 
list. This list consists of those vertex-edge pairs whose 
distance between them is less than + rj + 2a, where 
ri and rj are the sphero-radii. In each time step, only 
these vertex-edge pairs are involved in the contact force 
calculations. Overall, the Verlet list for neighborhood 
identification consists of a neighbor list with all pair of 
neighbor particles, and one contact list for each pair of 
neighbors. These lists require little memory storage, and 
they reduce the amount of calculations of contact forces 
to 0{N)^ which is of the same order as in simulations 
with spherical particles [11]. 

The Verlet list is calculated at the beginning of the sim- 
ulation, and it is updated when the following condition 
is satisfied: 



max {Axi 

l<i<N 



RiAOi} > a. 



(17) 



where fi is the center of mass of the particle i and R{V^ E) 
is the point of application of the force and . This is the 



Axi and A^^ are the maximal displacement and rotation 
of the particle after the last neighbor list update. Ri the 
maximal distance from the points on the particle to its 
center of mass. After each update Axi and AOi are set 
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to zero. The update condition is checked in each time 
step. 

Increasing the value of a makes updating of the hst less 
frequent, but increases its size, and hence the number of 
force calculations. Therefore, the parameter a must be 
chosen by chosing a reasonable balance between the time 
expended by the force calculations and the updating of 
neighbor lists. 

IV. SIMULATIONS 

In order to check the validity of this method we created 
a shear band simulation for our spheropolygons with pe- 
riodic boundary conditions. The polygons are confined 
by two horizontal plates. The lower plate is fixed. The 
upper plate has a fixed load and a constant horizontal ve- 
locity. First he verify energy balance by calculating the 
work done by the external forces, the dissipated energy, 
and the internal energy of the systems. Then we per- 
form a series of simulations to analize the computational 
efficiency of the method. 




FIG. 9: Shear band simulator with regular pentagons (above) 
and voronoi polygons (below) 



A. Energy Balance Test 

In this section we present numerical results from sim- 
ulations with dissipative granular materials. We check 
whether the system complies with the energy balance as 
established by the first law of the thermodynamics (work 



done by the system = energy dissipated + change of in- 
ernal energy of the system). 

Simulations are performed using two different particle 
geometries: regular spheropolyhedra and particles gen- 
erated with Voronoi-Minkowski diagrams, see Fig. [9] . 
These systems are simulated by replacing in Eq. [13] the 
following repulsive, frictional force 



F{V,E)=knSnN^kt6tf, (18) 



where N = ^ is the normal unit vector. The tan- 

gential vector T is taken perpendicular to N. The over- 
lapping length Sn is defined in Eq. [161 The elastic dis- 
placement St accounts frictional forces, and it must sat- 
isfied the sliding condition \Ft\ < jiFn^ where /i is the 
coefficient of friction [19*]. The parameters of the simula- 
tions are the normal stiffness k = 10000, the tangential 
stiffness kt = O.l/c^, friction coefficient /i = 0.5, density 
(7 = 1, and the time step At = 0.00001. Viscosity forces 
proportional to relative velocity of the contacts are also 
included to allows relaxation of the system. 

The evolution of and the orientation (pi of the par- 
ticle is obtained by solving the equations of motion: 



rriin ^Fji - ruigy, ^(pi ^r^i- (19) 

3 3 



Here rrii and li are the mass and moment of inertia of 
the particle. The sum is over all particles interacting 
with this particle; ^ = 10 is the gravity; and y is the 
unit vector along the vertical direction. The equations of 
motion of the system are numerically solved using a four 
order predictor-corrector algorithm . 

For both cases the total mechanical energy is obtained 
as well as the work done by the external forces. We show 
the relative difference for both cases (i.e. the difference 
between the mechanical energy plus de dissipated energy 
and the work done by external forces relative to the ex- 
ternal work) in Fig. [TOl 
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FIG. 10: Relative difference between the total mechanical 
energy and the dissipated one with the work by external forces 
for regular pentagons (above) and voronoi polygons (below) 



As can be seen the difference is less than 1% for both 
regular pentagons and voronoi polygons. The simulation 
is ready for problems that require energy measurement 
like the heating of internal water in the void spaces of 
the shear band or the increase of the rock temperature 
due to the frictional force. This is an important advan- 
tage of our model since the repulsive force as calculated 
provides a simple potential energy expression which was 
not available in previous models [ll|,ll2|. 



B. Performance Test 

One of the main optimization methods implemented 
in our code is the Verlet list. Certainly this method 
greatly reduce the amount of CPU time expended on 
the force calculations. The efficiency of the algorithm is 
meassure using the so-called Cundall numher[21\. This 
number is the amount of particle time steps executed 
by the processor in one second, which is calculated as 
c = NtN/Tcpui where Nt is the number of time steps, 
N is the number of particles and Tcpu is the CPU time 
of the simulation. The inverse of the Cundall number 1/c 
will be called here normalized CPU time Previous sim- 
ulations [21] shown that the Cundall number, although 
dependent on the processor, does not depend much on the 
number of particle if N > 100. Therefore for the simula- 
tion presented in this paper we chose samples with 102 
particles. 

Our first test for this optimization method was to mea- 
sure the CPU time of a shear simulation for grains of 
different number of sides. We start with triangles and 
end with heptagons, the results are shown in Fig. [TTl 
Initially with an a = 0.9 the Verlet list is not frequently 
updated and hence the optimization is for all purposes 
turned off. The force is then the predominant load on 
the CPU time and its calculation grows like a quadratic 
law with the numbers of sides. 
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FIG. 11: Inverse of the Cundall Number versus number of 
vertices of the constituting element. 



This change with the introduction of a smaller a value. 
The optimization reduce the force load on the total calcu- 
lation time. The power law dependence observed in the 
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previous case is changed (exponent is equal to 0.3533), 
and there is no apparent difference between the time 
expended in the calculation for hexagons and for hep- 
tagons. With the optimization the time no longer grows 
in a quadratic form. 

The next test was to measure the specific time that the 
simulation assign in both calculating forces and updating 
the Verlet list. The results for disks are show in Fig. [T2l 
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FIG. 12: The distinct times taken for a shear band simula- 
tion with disks elements with the profiling utility of the gcc 
compiler versus the a parameter 



The total CPU time is mainly split in force calculation 
and updating the Verlet lists. The time spent in each 
one of these features depends on a. Now for a small a 
the Verlet list are more frequently updated, and actually 
for a ^ the updating frequency will diverge to oo also. 
In contrast when a is large, the simulation no longer 
updates the Verlet list and therefore the time spent in 
this action goes to zero. Hence an invert proportional 
relation should be adequate to model this part of the 
total time. Now the time spent in the force calculation 
grows with a. We propose a quadratic law since the 
number of grains in a Verlet distance to the particle will 
be proportional to its area. 

With this idea me made the final test to our method: 
The construction of a shear band simulation (see Fig. 
[9]) in which the load is moving to the right side with a 
constant velocity. We check the CPU time vs the velocity 
of the load for a range of a values going from 0.1 to 0.9 
and with three different load velocities (0.1, 0.4,0.7 and 
1.). The results are shown in Fig. [131 

The results where fitted with this analytical model re- 
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FIG. 13: Obtained cpu time for the shear band test with 
different shear rates. The fitted curves are obtained with the 
model proposed on Eq. [20] 



lating the CPU time with the verlet length a: 



Tc 



PU 



'-f 



A 



(20) 



Where T^c is the time spent on updating the verlet 
list, Tf is the time for calculating the forces, A, B and 
P are de fitting parameters. As we can see the optimum 
value value for a change with the shear rate as well as the 
total CPU time for a given a. For the lowest shear rate 
we have the best a value, which means that the Verlet 
list method accomplishes a greater optimization if the 
particles move slowly. 



V. CONCLUDING REMARKS 

The main motivation of this paper was to show that 
the concepts of computational geometry (Voronoi dia- 
grams, minkowski operator and distance formulas) are 
useful to model systems composed with complex shaped 
particles. These concepts are applied to simulate systems 
of spheropolygons interacting via dissipative forces. We 
guarantee energy balance, which were not available in 
previous models. The classical concept of Verlet list was 
extended for spheropolygons. Depending on the charac- 
teristic velocity of the system there is an optimal value 
of the Verlet distance that minimize CPU time. 

The most promising aspect of this work is its natural 
extension to 3D. There are several available subroutines 
to perform Voronoi Tesselation in 3D, which will allow to 



9 



generate the Voronoi-Minkowski diagrams once the sub- 
routine of erosion is implemented. The analytical calcula- 
tion of the mass properties in 3D spheropolyhedra can be 
obtained by decomposing them in polyhedra and spher- 
ical segments. The contact force between spheropolyhe- 
dra can be calculated as the sum of all vertex-face and 
edge-edge contacts. The equation of motion for the ro- 
tational degrees of freedom are not too trivial like in 2D, 
but the quaternion formalism [22] can be used to a cor- 



rect description of all the degrees of freedom. 
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